close all
clear all
clc
format long
global a; global b; global c;
abc=[1 1 1;1 0 100; -1 0 100; -10 10 0; 1 10000 0; 1 10000 250; 10000 1 250];
tolerances=[1e-3 1e-8];
figure1=figure('Name','relatieve tolerantie 10e-3');
figure2=figure('Name','relatieve tolerantie 10e-8');
figures=[figure1 figure2];

feval = [];
avgstep = [];
for i=1:2,
    
    Tol=tolerances(i);
    options = odeset('RelTol',Tol,'Stats','on');
    for j=1:7,
        
        %% Instellen van de parameters
        a=abc(j,1);
        b=abc(j,2);
        c=abc(j,3);
        N=50;
        tspan=[0:1/N:1]';
        
        %% Beginvoorwaarden
        
        %   y1(0)=2
        %   y2(0)=c-a
        %   y3(0)=a^2-2bc
        
        y0=[2 c-a a^2-2*b*c];
        
        %% Opstellen van het stelsel
        %
        %
        % met y(t) = [y1(t) y2(t) y3(t)]
        % en y1'(t)=y2(t)
        %    y2'(t)=y3(t)
        %    y3'(t)=a(b^2+c^2)-(2b+a)y3(t)-(b^2+c^2+2ab)y2(t)-a(b^2+c^2)y1(t)
        %
        % matrix A bepalen uit y'(x)=Ay(x)+f(t);
        
        A=[0 1 0; 0 0 1; -a*(b^2+c^2) -(b^2+c^2+2*a*b) -(2*b+a)];
        D=eig(A);
        k=max(abs(real(eig(A))))/min(abs(real(eig(A))));
        display(['a = ' num2str(a) ', b = ' num2str(b) ', c= ' num2str(c) ', k= ' num2str(k) ', Tol= ' num2str(Tol)]); display(' ');
        eigenwaarden = eig(A)
        %% Oplossen van het stelsel: ode45 (Runge-Kutta)
        
        t=cputime;
        [tout1 yout1 stat]=ode45(@odefun,tspan,y0,options);
        stat
        feval = [feval;stat(3)];
        avgstep = [avgstep;1/stat(1)];
        display(['cputime = ' num2str(cputime-t)]); 
        fprintf('average step length %f \n',1/stat(1));display(' ')
        %% Oplossen van het stelsel: ode113 (VSVO-AB)
        
        t=cputime;
        [tout2 yout2 stat]=ode113(@odefun,tspan,y0,options);
        stat
        feval = [feval;stat(3)];
        avgstep = [avgstep;1/stat(1)];
        display(['cputime = ' num2str(cputime-t)]); 
        fprintf('average step length %f \n',1/stat(1));display(' ');
        
        %% Oplossen van het stelsel: ode15s (NDF)
        
        t=cputime;
        [tout3 yout3 stat]=ode15s(@odefun,tspan,y0,options);
        stat
        feval = [feval;stat(3)];
        avgstep = [avgstep;1/stat(1)];
        display(['cputime = ' num2str(cputime-t)]); 
        fprintf('average step length %f \n',1/stat(1));display(' ');
        
        display('-----------------------'); display(' ');
        
        %% Oplossen van het stelsel: exacte oplossing
        
        texact=tspan;
        yexact=1+exp(-b*tspan).*sin(c*tspan)+exp(-a*tspan);
        
        figure(figures(i));
        subplot(3,3,j);
        hold on
        plot(texact,-(log(abs(yout1(:,1)-yexact)./yexact)),'r','Linewidth',2)
        plot(texact,-(log(abs(yout2(:,1)-yexact)./yexact)),'b','Linewidth',2)
        plot(texact,-(log(abs(yout3(:,1)-yexact)./yexact)),'g','Linewidth',2)
        
        title(['a = ' num2str(a) ', b = ' num2str(b) ', c= ' num2str(c) ', k= ' num2str(k)]);
        
        legend('ode45','ode113','ode15s')
        hold off
        
        figure();
        hold on
        plot(texact,-(log(abs(yout1(:,1)-yexact)./yexact)),'r','Linewidth',2)
        plot(texact,-(log(abs(yout2(:,1)-yexact)./yexact)),'b','Linewidth',2)
        plot(texact,-(log(abs(yout3(:,1)-yexact)./yexact)),'g','Linewidth',2)
        xlabel('x');
        ylabel('Relatieve fout');
        legend('ode45','ode113','ode15s')
        % matlab2tikz(sprintf('RelErrTol%dSet%d.tikz', i, j), 'height', '\figureheight', 'width', '\figurewidth' );
        close gcf
    end
end

feval
avgstep
